Assignment 2: Multi-View Geometry

Comment: most of the written problems (except the first one) are designed to help with the coding part (structure-from-motion). Thus, they should be solved first. Your coding part require slicing, broadcasting, and other standard operations with numpy arrays. This could be tricky if you are new to numpy. I recommend working on small test cases. For debugging, print arrays before and after slicing (etc.) to verify that the result is correct.

Problem 1 (Frobenius norm)

Compute Frobenius norm of $n\times n$ matrix $xx^\top$ assuming that (Euclidean) norm $\|x\| = \sqrt{x^\top x}$ of vector $x\in R^n$ is given. That is, express $\|xx^\top\|_F$ in terms of $\|x\|$. The answer should not depend on $n$. Comment: this "random math" excerise is only meant to encourage you to look up the definition of Frobenius norm; the proof is only a couple of lines of simple algebra.

Solution:

$$ xx^T = \begin{bmatrix} x_1^2 & x_1x_2 & \cdots & x_1x_n \\ x_2x_1 & x_2^2 & \cdots & x_2x_n \\ \vdots & \vdots & \ddots & \vdots \\ x_nx_1 & x_nx_2 & \cdots & x_n^2 \end{bmatrix} $$$$ \|xx^T\|_F = \sqrt{trace((xx^T)^T(xx^T))} \\ = \sqrt{trace((x^Tx)(xx^T))} \\ = \sqrt{trace(\|x\|^2(xx^T))} \\ = \|x\|\sqrt{trace(xx^T)} \\ = \|x\|\sqrt{\|x\|^2} \\ = \|x\|^2 $$

Problem 2

Assuming a $calibrated$ camera (that is, $K=I$) and its two views corresponding to projection matrices $P_1=[I|0]$ and $P_2=[R|T]$ w.r.t. some world coordinate system, show formulas for coordinates of the following 3D points (in the same world coordinate system):

$$ \begin{bmatrix} x \\ y \\ z \end{bmatrix} = P \begin{bmatrix} X \\ Y \\ Z \\ 1 \end{bmatrix} \\ \begin{bmatrix} wu \\ wv \\ w \end{bmatrix} = KP \begin{bmatrix} X \\ Y \\ Z \\ 1 \end{bmatrix} $$

(a) optical center for the first view: $C_1= (0,0,0)$

(b) image center for the first view: $Q_1= (0, 0, 1)$

(c) optical center for the second view: $C_2=(RC_1 + T) = T$

(d) image center for the second view: $Q_2= (RQ_1 + T)$

Problem 3

Using the same set up as in problem 2, show formulas for normalized coordinates of the following image points:

(a) epipole in the first camera image: $e_1=R_1C_2 + T_1 = C_2 = T$

(b) epipole in the second camera image: $e_2=RC_1 + T = T$

Problem 4 (homogeneous and non-homogeneous line representations)

Lines in 2D images can be represented "homogeneously" as 3-vectors $l=[l_1,l_2,l_3]^T$ that give equation $l^T x=0$ for homogeneous points $x=[x_1,x_2,x_3]^T \;\in {\cal P}^2$ forming a line. Given $l$, what are the values of scalar parameters $a$, $b$ in the line equation $v=au+b$ for the same 2D points based on their regular (nonhomogeneous) representation $(u,v)=(\frac{x_1}{x_3},\frac{x_2}{x_3})$ in ${\cal R}^2$?

$$ l^Tx = 0 \\ x_1l_1+x_2l_2+x_3l_3 = 0 \\ \frac{x_1}{x_3}\frac{l_1}{l_2} + \frac{x_2}{x_3} + \frac{l_3}{l_2} = 0 \\ \frac{x_2}{x_3} = (\frac{-l_1}{l_2})\frac{x_1}{x_2} + \frac{-l_3}{l_2} \\ v = au + b $$

$a=\frac{-l_1}{l_2}$

$b=\frac{-l_3}{l_2}$

Problem 5 (epipolar lines in normalized and non-normalized images)

Given a matrix of intrinsic camera parameters $K$ and essential matrix $E$ between two views (A) and (B) such that $x_A^T E x_B=0$ for any corresponding points, write expressions for the following:

(a) given homogeneous normalized point $x^{n}_B$ in image B, specify 3-vector $l_A^n$ describing the corresponding epipolar line of normalized points in image A:

$l_A^n=E^Tx^n_B$

(b) given homogeneous normalized point $x^{n}_A$ in image A, specify 3-vector $l_B^n$ describing the corresponding epipolar line of normalized points in image B:

$l_B^n=Ex^n_A$

(c) assuming line (3-vector) $l^n$ of normalized image points, what is a 3-vector representation $l$ for the line formed by the corresponding points on the real (unnormalized) camera image:

$l=Kl^n$

Problem 6 (least squares for triangulation)

Describe your approach to triangulating two matched feature points $x_a=[u_a,v_a,1]^T$ and $x_b=[u_b,v_b,1]^T$ in two views with given projection matrices $P_a$ and $P_b$. You should find 3D point $X=[X_1,X_2,X_3,1]^T$ and two scalars $w_a,w_b$ such that $P_a X\approx w_a x_a$ and $P_b X\approx w_b x_b$. Be specific as you will need this for your programming part below. Use notation $M[i]$ to denote the $i$-th row vector of matrix $M$.

You should use the first approach described for homograpy estimation in topic 6. In particualr, you can formulate the problem as $AX\approx 0$, define elements of $4x4$ matrix $A$, convert the problem to an overdetermined system of 4 linear equations $A_{1:3}[X_1,X_2,X_3]^T \approx - A_{4}$, and specify its solution minimizing the sum of squared errors.

Can you characterize geometrically the case when your solution satisfies $A_{1:3}[X_1,X_2,X_3]^T = - A_{4}$ exactly?

Solution:

Given that $w_ax_a \approx P_aX$, we can use a cross product to remove the scalar term s.t. $$ \begin{bmatrix} u_a \\ v_a \\ 1 \end{bmatrix} \times P_aX \\ = \begin{bmatrix} u_a \\ v_a \\ 1 \end{bmatrix} \times \begin{bmatrix} p_1^TX \\ p_2^TX \\ p_3^TX \end{bmatrix} \\ = \begin{bmatrix} v_ap_3^TX - p_2^TX \\ p_1^TX - xp_3^TX \\ xp_2^TX - yp_1^TX \end{bmatrix} \\ = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} $$

We can construct A using above for each point, resulting in: $$ A = \begin{bmatrix} v_ap_{a3}^TX - p_{a2}^TX \\ p_{a1}^TX - u_ap_{a3}^TX \\ v_ap_{b3}^TX - p_{b2}^TX \\ p_{b1}^TX - u_ap_{b3}^TX \end{bmatrix} $$

Geometrically, solving for A_{1:3} = -A_4 means that the solutions are on the epipolar lines.

Probelm 7 (the programming part)

Structure from Motion

NOTE: Steps 0-3 and 10 are given, other steps are to be implemented.

Step 0: Loading two camera views and camera's intrinsic matrix $K$

%matplotlib inline

import numpy as np
import numpy.linalg as la
import matplotlib
import matplotlib.image as image
import matplotlib.pyplot as plt
from skimage.feature import (corner_harris, corner_peaks, plot_matches, BRIEF, match_descriptors)
from skimage.transform import warp, ProjectiveTransform, EssentialMatrixTransform, FundamentalMatrixTransform
from skimage.color import rgb2gray
from skimage.measure import ransac

# Indicate (E) inlier matches in image 1 and image 2
# loading two images (two camera views) and the corresponding matrix K (intrinsic parameters)
imL = image.imread("images/kronan1.jpg")
imR = image.imread("images/kronan2.jpg")
imLgray = rgb2gray(imL)
imRgray = rgb2gray(imR)

K = 1.0e+03 * np.array([[2.3940, -0.0000,    0.9324],
                        [     0,  2.3981,    0.6283],
                        [     0,       0,    0.0010]])


plt.figure(0,figsize = (10, 4))
ax81 = plt.subplot(121)
plt.imshow(imL)
ax82 = plt.subplot(122)
plt.imshow(imR)
plt.show()

Step 1: Feature detection (e.g. corners)

# NOTE: corner_peaks and many other feature extraction functions return point coordinates as (y,x), that is (rows,cols)
keypointsL = corner_peaks(corner_harris(imLgray), threshold_rel=0.001, min_distance=15)
keypointsR = corner_peaks(corner_harris(imRgray), threshold_rel=0.001, min_distance=15)


print ('the number of features in images 1 and 2 are {:5d} and {:5d}'.format(keypointsL.shape[0],keypointsR.shape[0]))

fig = plt.figure(1,figsize = (10, 4))
axA = plt.subplot(111)
plt.gray()
matchesLR = np.empty((0,2))
plot_matches(axA, imL, imR, keypointsL, keypointsR, matchesLR)
axA.axis('off')

plt.show()

the number of features in images 1 and 2 are  1576 and  1661

Step 2: Feature matching (e.g. BRIEF descriptor, a variant of SURF, SIFT, etc)

extractor = BRIEF()

extractor.extract(imLgray, keypointsL)
keypointsL = keypointsL[extractor.mask]         
descriptorsL = extractor.descriptors

extractor.extract(imRgray, keypointsR)
keypointsR = keypointsR[extractor.mask]
descriptorsR = extractor.descriptors

matchesLR = match_descriptors(descriptorsL, descriptorsR, cross_check=True)

print ('the number of matches is {:2d}'.format(matchesLR.shape[0]))

fig = plt.figure(2,figsize = (10, 4))
axA = plt.subplot(111)
axA.set_title("matches")
plt.gray()
plot_matches(axA, imL, imR, keypointsL, keypointsR, matchesLR) #, matches_color = 'r')
axA.axis('off')

plt.show()

the number of matches is 964

Step 3: Fundamental Matrix estimation using RANSAC

ptsL1 = []
ptsR1 = []
for i in matchesLR:
    ptsL1.append(keypointsL[i[0]])
    ptsR1.append(keypointsR[i[1]])
ptsL1 = np.array(ptsL1)
ptsR1 = np.array(ptsR1)

# swapping columns using advanced indexing https://docs.scipy.org/doc/numpy/reference/arrays.indexing.html#advanced-indexing
# This changes point coordinates from (y,x) in ptsL1/ptsR1 to (x,y) in ptsL/ptsR
ptsL = ptsL1[:,[1, 0]]
ptsR = ptsR1[:,[1, 0]]

# robustly estimate fundamental matrix using RANSAC
F_trans, F_inliers = ransac((ptsL, ptsR), FundamentalMatrixTransform, min_samples=8, residual_threshold=0.1, max_trials=1500)
print ('the number of inliers is {:2d}'.format(np.sum(F_inliers)))

ind = np.ogrid[:ptsL.shape[0]]
FmatchesRansac = np.column_stack((ind[F_inliers],ind[F_inliers]))

fig = plt.figure(3,figsize = (10, 4))
axA = plt.subplot(111)
axA.set_title("inlier matches")
plt.gray()
# NOTE: function "plot matches" expects that keypoint coordinates are given as (y,x), that is (row, col)
plot_matches(axA, imL, imR, ptsL1, ptsR1, FmatchesRansac) #, matches_color = 'r')
axA.axis('off')
plt.show()

the number of inliers is 178

singular values for F

F = F_trans.params
Uf,Sf,Vf = la.svd(F, full_matrices=False)
print (Sf)

[7.73816619e-02 6.00151466e-05 2.52310535e-20]

Step 4: Epipolar lines from F

# Randomly select 10 matches (paris of features in two images) from the set of inliers for F
ind_sample = np.random.choice(ind[F_inliers], 10, replace = False)

# Indicate these matching features in image 1 and image 2
plt.figure(4,figsize = (10, 4))
ax41 = plt.subplot(121)
plt.imshow(imL)
plt.plot(ptsL[ind_sample, 0], ptsL[ind_sample, 1], 'ob')
ax42 = plt.subplot(122)
plt.imshow(imR)
plt.plot(ptsR[ind_sample, 0], ptsR[ind_sample, 1], 'ob')

# generate epipolar line equations in image 2 (homoheneous 3-vectors l2 representing lines l2 x  = 0) 
# a. create an array of points sampled in images 1 and 2  
# b. create an array of homogeneous points sampled in images 1 and 2 
# c. create an array of the corresponding epipolar lines in images 1 and 2 
ones = np.ones(ptsL[ind_sample].shape[0])
pts_L_h = np.column_stack((ptsL[ind_sample, 0], ptsL[ind_sample, 1], ones))
pts_R_h = np.column_stack((ptsR[ind_sample, 0], ptsR[ind_sample, 1], ones))

# l_L = np.einsum('ij,kj->ki', F, pts_L_h)
# l_R = np.einsum('ij,kj->ki', F, pts_R_h)

for i in range(10):
    l_R = np.matmul(F, pts_L_h[i])
    l_L = np.matmul(F.T, pts_R_h[i])

    x = np.linspace(0, imL.shape[1])
    y_R = (-l_R[0]/l_R[1])*x + (-l_R[2]/l_R[1])
    y_L = (-l_L[0]/l_L[1])*x + (-l_L[2]/l_L[1])

    ax41.plot(x, y_L, 'k')
    ax42.plot(x, y_R, 'k')
    
# for each feature (in both images) draw a correspoindiung epipolar line in the other image
# see Assignment 1 (line fitting part 1) for inspiration on how to visualize lines
# use ax41.plot and ax42.plot 


plt.show()

Step 5: Camera Normalization and Essential Matrix estimation using RANSAC

# normalization of points in two images using K (intrinsic parameters) e.g. in the following three steps
# a. convert original points to homogeneous 3-vectors (append "1" as a 3rd coordinate using np.append function)
# b. transform the point by applying the inverse of K
# c. convert homogeneous 3-vectors to 2-vectors (in R2)
n_ptsL = np.column_stack((ptsL[:, 0], ptsL[:, 1], np.ones(ptsL.shape[0])))
n_ptsR = np.column_stack((ptsR[:, 0], ptsR[:, 1], np.ones(ptsR.shape[0])))


K_inv = np.linalg.inv(K)
tr_ptsL_0 = np.matmul(K_inv, n_ptsL.T).T
tr_ptsR_0 = np.matmul(K_inv, n_ptsR.T).T


tr_ptsL = np.einsum('ij,kj->ki', K_inv, n_ptsL)
tr_ptsR = np.einsum('ij,kj->ki', K_inv, n_ptsR)

print(tr_ptsL - tr_ptsL_0)

# print(tr_ptsL[:, 2].shape)

n_ptsL = np.divide(tr_ptsL[:, :2], tr_ptsL[:, 2][:, np.newaxis])
n_ptsR = np.divide(tr_ptsL[:, :2], tr_ptsL[:, 2][:, np.newaxis])

# print(n_ptsL)

# robustly estimate essential matrix using normalized points and RANSAC
E_trans, E_inliers = ransac((n_ptsL, n_ptsR), EssentialMatrixTransform, min_samples=8, residual_threshold=0.0005, max_trials=5000)
num_inliers = np.sum(E_inliers)
print ('the number of inliers is {:2d}'.format(num_inliers))

ind = np.ogrid[:n_ptsL.shape[0]]
EmatchesRansac = np.column_stack((ind[E_inliers],ind[E_inliers]))

fig = plt.figure(5,figsize = (10, 4))
axA = plt.subplot(111)
axA.set_title("inlier matches")
plt.gray()
# NOTE: function "plot matches" expects that keypoint coordinates are given as (y,x), that is (row, col)
plot_matches(axA, imL, imR, ptsL1, ptsR1, EmatchesRansac) #, matches_color = 'r')
axA.axis('off')
plt.show()

[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]
 ...
 [0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
the number of inliers is 964

singular values for E

Hint: function $svd$ from $linalg$ returns transpose $V^T$, not $V$.

E = E_trans.params
Ue,Se,Ve = la.svd(E)
print (Se)

[2.53153836e+01 2.53153836e+01 6.84258394e-18]
(3, 3)

Step 6: Epipolar Lines from E

# Randomly select 10 matches (paris of features in two images) from the set of inliers for E
ind_sample = np.random.choice(ind[E_inliers], 10, replace = False)

# Indicate these matching features in image 1 and image 2
plt.figure(6,figsize = (10, 4))
ax61 = plt.subplot(121)
plt.imshow(imL)
plt.plot(ptsL[ind_sample, 0], ptsL[ind_sample, 1], 'ob')
ax62 = plt.subplot(122)
plt.imshow(imR)
plt.plot(ptsR[ind_sample, 0], ptsR[ind_sample, 1], 'ob')

# generate epipolar line equations in image 2 (homoheneous 3-vectors l2 representing lines l2 x  = 0) 
# a. create an array of normalized points sampled in image 1 
# b. create an array of homogeneous normalized points sampled in image 1 
# c. create an array of the corresponding (uncalibrated) epipolar lines in image 2 
n_ptsL_sample = np.column_stack((n_ptsL[ind_sample, 0], n_ptsL[ind_sample, 1], np.ones(n_ptsL[ind_sample].shape[0])))
n_ptsR_sample = np.column_stack((n_ptsR[ind_sample, 0], n_ptsR[ind_sample, 1], np.ones(n_ptsR[ind_sample].shape[0])))

x = np.arange(0, imL.shape[1], step=5)

for i in range(10):
    l_R = np.matmul(E, n_ptsL_sample[i])
    l_L = np.matmul(E.T, n_ptsR_sample[i])
    print(l_L)

    y_R = (-l_R[0]/l_R[1])*x + (-l_R[2]/l_R[1])
    y_L = (-l_L[0]/l_L[1])*x + (-l_L[2]/l_L[1])

    m_L = ((y_L < imLgray.shape[0]) & (y_L >= 0)).any()
    m_R = ((y_R < imLgray.shape[0]) & (y_R >= 0)).any()

    print(y_R[m_R])

    ax61.plot(x, y_L, 'k')
    ax62.plot(x, y_R, 'k')

# for each feature (in both images) draw a correspoindiung epipolar line in the other image
# use ax61.plot and ax62.plot 

plt.show()

[-0.46242028  3.71029245 -0.11726585]
[[3.16055545e-02 6.54764357e-01 1.27792316e+00 1.90108196e+00
  2.52424077e+00 3.14739957e+00 3.77055837e+00 4.39371718e+00
  5.01687598e+00 5.64003478e+00 6.26319358e+00 6.88635239e+00
  7.50951119e+00 8.13266999e+00 8.75582880e+00 9.37898760e+00
  1.00021464e+01 1.06253052e+01 1.12484640e+01 1.18716228e+01
  1.24947816e+01 1.31179404e+01 1.37410992e+01 1.43642580e+01
  1.49874168e+01 1.56105756e+01 1.62337344e+01 1.68568932e+01
  1.74800520e+01 1.81032108e+01 1.87263696e+01 1.93495284e+01
  1.99726872e+01 2.05958461e+01 2.12190049e+01 2.18421637e+01
  2.24653225e+01 2.30884813e+01 2.37116401e+01 2.43347989e+01
  2.49579577e+01 2.55811165e+01 2.62042753e+01 2.68274341e+01
  2.74505929e+01 2.80737517e+01 2.86969105e+01 2.93200693e+01
  2.99432281e+01 3.05663869e+01 3.11895457e+01 3.18127045e+01
  3.24358633e+01 3.30590221e+01 3.36821809e+01 3.43053397e+01
  3.49284985e+01 3.55516573e+01 3.61748161e+01 3.67979749e+01
  3.74211337e+01 3.80442925e+01 3.86674513e+01 3.92906101e+01
  3.99137689e+01 4.05369277e+01 4.11600866e+01 4.17832454e+01
  4.24064042e+01 4.30295630e+01 4.36527218e+01 4.42758806e+01
  4.48990394e+01 4.55221982e+01 4.61453570e+01 4.67685158e+01
  4.73916746e+01 4.80148334e+01 4.86379922e+01 4.92611510e+01
  4.98843098e+01 5.05074686e+01 5.11306274e+01 5.17537862e+01
  5.23769450e+01 5.30001038e+01 5.36232626e+01 5.42464214e+01
  5.48695802e+01 5.54927390e+01 5.61158978e+01 5.67390566e+01
  5.73622154e+01 5.79853742e+01 5.86085330e+01 5.92316918e+01
  5.98548506e+01 6.04780094e+01 6.11011682e+01 6.17243270e+01
  6.23474859e+01 6.29706447e+01 6.35938035e+01 6.42169623e+01
  6.48401211e+01 6.54632799e+01 6.60864387e+01 6.67095975e+01
  6.73327563e+01 6.79559151e+01 6.85790739e+01 6.92022327e+01
  6.98253915e+01 7.04485503e+01 7.10717091e+01 7.16948679e+01
  7.23180267e+01 7.29411855e+01 7.35643443e+01 7.41875031e+01
  7.48106619e+01 7.54338207e+01 7.60569795e+01 7.66801383e+01
  7.73032971e+01 7.79264559e+01 7.85496147e+01 7.91727735e+01
  7.97959323e+01 8.04190911e+01 8.10422499e+01 8.16654087e+01
  8.22885675e+01 8.29117263e+01 8.35348852e+01 8.41580440e+01
  8.47812028e+01 8.54043616e+01 8.60275204e+01 8.66506792e+01
  8.72738380e+01 8.78969968e+01 8.85201556e+01 8.91433144e+01
  8.97664732e+01 9.03896320e+01 9.10127908e+01 9.16359496e+01
  9.22591084e+01 9.28822672e+01 9.35054260e+01 9.41285848e+01
  9.47517436e+01 9.53749024e+01 9.59980612e+01 9.66212200e+01
  9.72443788e+01 9.78675376e+01 9.84906964e+01 9.91138552e+01
  9.97370140e+01 1.00360173e+02 1.00983332e+02 1.01606490e+02
  1.02229649e+02 1.02852808e+02 1.03475967e+02 1.04099126e+02
  1.04722284e+02 1.05345443e+02 1.05968602e+02 1.06591761e+02
  1.07214920e+02 1.07838078e+02 1.08461237e+02 1.09084396e+02
  1.09707555e+02 1.10330714e+02 1.10953872e+02 1.11577031e+02
  1.12200190e+02 1.12823349e+02 1.13446508e+02 1.14069666e+02
  1.14692825e+02 1.15315984e+02 1.15939143e+02 1.16562302e+02
  1.17185461e+02 1.17808619e+02 1.18431778e+02 1.19054937e+02
  1.19678096e+02 1.20301255e+02 1.20924413e+02 1.21547572e+02
  1.22170731e+02 1.22793890e+02 1.23417049e+02 1.24040207e+02
  1.24663366e+02 1.25286525e+02 1.25909684e+02 1.26532843e+02
  1.27156001e+02 1.27779160e+02 1.28402319e+02 1.29025478e+02
  1.29648637e+02 1.30271795e+02 1.30894954e+02 1.31518113e+02
  1.32141272e+02 1.32764431e+02 1.33387589e+02 1.34010748e+02
  1.34633907e+02 1.35257066e+02 1.35880225e+02 1.36503383e+02
  1.37126542e+02 1.37749701e+02 1.38372860e+02 1.38996019e+02
  1.39619177e+02 1.40242336e+02 1.40865495e+02 1.41488654e+02
  1.42111813e+02 1.42734971e+02 1.43358130e+02 1.43981289e+02
  1.44604448e+02 1.45227607e+02 1.45850765e+02 1.46473924e+02
  1.47097083e+02 1.47720242e+02 1.48343401e+02 1.48966559e+02
  1.49589718e+02 1.50212877e+02 1.50836036e+02 1.51459195e+02
  1.52082353e+02 1.52705512e+02 1.53328671e+02 1.53951830e+02
  1.54574989e+02 1.55198147e+02 1.55821306e+02 1.56444465e+02
  1.57067624e+02 1.57690783e+02 1.58313942e+02 1.58937100e+02
  1.59560259e+02 1.60183418e+02 1.60806577e+02 1.61429736e+02
  1.62052894e+02 1.62676053e+02 1.63299212e+02 1.63922371e+02
  1.64545530e+02 1.65168688e+02 1.65791847e+02 1.66415006e+02
  1.67038165e+02 1.67661324e+02 1.68284482e+02 1.68907641e+02
  1.69530800e+02 1.70153959e+02 1.70777118e+02 1.71400276e+02
  1.72023435e+02 1.72646594e+02 1.73269753e+02 1.73892912e+02
  1.74516070e+02 1.75139229e+02 1.75762388e+02 1.76385547e+02
  1.77008706e+02 1.77631864e+02 1.78255023e+02 1.78878182e+02
  1.79501341e+02 1.80124500e+02 1.80747658e+02 1.81370817e+02
  1.81993976e+02 1.82617135e+02 1.83240294e+02 1.83863452e+02
  1.84486611e+02 1.85109770e+02 1.85732929e+02 1.86356088e+02
  1.86979246e+02 1.87602405e+02 1.88225564e+02 1.88848723e+02
  1.89471882e+02 1.90095040e+02 1.90718199e+02 1.91341358e+02
  1.91964517e+02 1.92587676e+02 1.93210834e+02 1.93833993e+02
  1.94457152e+02 1.95080311e+02 1.95703470e+02 1.96326628e+02
  1.96949787e+02 1.97572946e+02 1.98196105e+02 1.98819264e+02
  1.99442423e+02 2.00065581e+02 2.00688740e+02 2.01311899e+02
  2.01935058e+02 2.02558217e+02 2.03181375e+02 2.03804534e+02
  2.04427693e+02 2.05050852e+02 2.05674011e+02 2.06297169e+02
  2.06920328e+02 2.07543487e+02 2.08166646e+02 2.08789805e+02
  2.09412963e+02 2.10036122e+02 2.10659281e+02 2.11282440e+02
  2.11905599e+02 2.12528757e+02 2.13151916e+02 2.13775075e+02
  2.14398234e+02 2.15021393e+02 2.15644551e+02 2.16267710e+02
  2.16890869e+02 2.17514028e+02 2.18137187e+02 2.18760345e+02
  2.19383504e+02 2.20006663e+02 2.20629822e+02 2.21252981e+02
  2.21876139e+02 2.22499298e+02 2.23122457e+02 2.23745616e+02
  2.24368775e+02 2.24991933e+02 2.25615092e+02 2.26238251e+02
  2.26861410e+02 2.27484569e+02 2.28107727e+02 2.28730886e+02
  2.29354045e+02 2.29977204e+02 2.30600363e+02 2.31223521e+02
  2.31846680e+02 2.32469839e+02 2.33092998e+02 2.33716157e+02
  2.34339315e+02 2.34962474e+02 2.35585633e+02 2.36208792e+02
  2.36831951e+02 2.37455109e+02 2.38078268e+02 2.38701427e+02
  2.39324586e+02 2.39947745e+02 2.40570904e+02 2.41194062e+02]]
[ 1.44669991  2.04091353 -0.01365405]
[[ 6.69016421e-03 -3.53755581e+00 -7.08180179e+00 -1.06260478e+01
  -1.41702937e+01 -1.77145397e+01 -2.12587857e+01 -2.48030317e+01
  -2.83472776e+01 -3.18915236e+01 -3.54357696e+01 -3.89800156e+01
  -4.25242615e+01 -4.60685075e+01 -4.96127535e+01 -5.31569995e+01
  -5.67012454e+01 -6.02454914e+01 -6.37897374e+01 -6.73339834e+01
  -7.08782293e+01 -7.44224753e+01 -7.79667213e+01 -8.15109673e+01
  -8.50552133e+01 -8.85994592e+01 -9.21437052e+01 -9.56879512e+01
  -9.92321972e+01 -1.02776443e+02 -1.06320689e+02 -1.09864935e+02
  -1.13409181e+02 -1.16953427e+02 -1.20497673e+02 -1.24041919e+02
  -1.27586165e+02 -1.31130411e+02 -1.34674657e+02 -1.38218903e+02
  -1.41763149e+02 -1.45307395e+02 -1.48851641e+02 -1.52395887e+02
  -1.55940133e+02 -1.59484379e+02 -1.63028625e+02 -1.66572871e+02
  -1.70117117e+02 -1.73661363e+02 -1.77205609e+02 -1.80749855e+02
  -1.84294101e+02 -1.87838347e+02 -1.91382593e+02 -1.94926838e+02
  -1.98471084e+02 -2.02015330e+02 -2.05559576e+02 -2.09103822e+02
  -2.12648068e+02 -2.16192314e+02 -2.19736560e+02 -2.23280806e+02
  -2.26825052e+02 -2.30369298e+02 -2.33913544e+02 -2.37457790e+02
  -2.41002036e+02 -2.44546282e+02 -2.48090528e+02 -2.51634774e+02
  -2.55179020e+02 -2.58723266e+02 -2.62267512e+02 -2.65811758e+02
  -2.69356004e+02 -2.72900250e+02 -2.76444496e+02 -2.79988742e+02
  -2.83532988e+02 -2.87077234e+02 -2.90621480e+02 -2.94165726e+02
  -2.97709972e+02 -3.01254218e+02 -3.04798464e+02 -3.08342710e+02
  -3.11886956e+02 -3.15431202e+02 -3.18975448e+02 -3.22519694e+02
  -3.26063940e+02 -3.29608186e+02 -3.33152432e+02 -3.36696678e+02
  -3.40240923e+02 -3.43785169e+02 -3.47329415e+02 -3.50873661e+02
  -3.54417907e+02 -3.57962153e+02 -3.61506399e+02 -3.65050645e+02
  -3.68594891e+02 -3.72139137e+02 -3.75683383e+02 -3.79227629e+02
  -3.82771875e+02 -3.86316121e+02 -3.89860367e+02 -3.93404613e+02
  -3.96948859e+02 -4.00493105e+02 -4.04037351e+02 -4.07581597e+02
  -4.11125843e+02 -4.14670089e+02 -4.18214335e+02 -4.21758581e+02
  -4.25302827e+02 -4.28847073e+02 -4.32391319e+02 -4.35935565e+02
  -4.39479811e+02 -4.43024057e+02 -4.46568303e+02 -4.50112549e+02
  -4.53656795e+02 -4.57201041e+02 -4.60745287e+02 -4.64289533e+02
  -4.67833779e+02 -4.71378025e+02 -4.74922271e+02 -4.78466517e+02
  -4.82010763e+02 -4.85555008e+02 -4.89099254e+02 -4.92643500e+02
  -4.96187746e+02 -4.99731992e+02 -5.03276238e+02 -5.06820484e+02
  -5.10364730e+02 -5.13908976e+02 -5.17453222e+02 -5.20997468e+02
  -5.24541714e+02 -5.28085960e+02 -5.31630206e+02 -5.35174452e+02
  -5.38718698e+02 -5.42262944e+02 -5.45807190e+02 -5.49351436e+02
  -5.52895682e+02 -5.56439928e+02 -5.59984174e+02 -5.63528420e+02
  -5.67072666e+02 -5.70616912e+02 -5.74161158e+02 -5.77705404e+02
  -5.81249650e+02 -5.84793896e+02 -5.88338142e+02 -5.91882388e+02
  -5.95426634e+02 -5.98970880e+02 -6.02515126e+02 -6.06059372e+02
  -6.09603618e+02 -6.13147864e+02 -6.16692110e+02 -6.20236356e+02
  -6.23780602e+02 -6.27324848e+02 -6.30869094e+02 -6.34413339e+02
  -6.37957585e+02 -6.41501831e+02 -6.45046077e+02 -6.48590323e+02
  -6.52134569e+02 -6.55678815e+02 -6.59223061e+02 -6.62767307e+02
  -6.66311553e+02 -6.69855799e+02 -6.73400045e+02 -6.76944291e+02
  -6.80488537e+02 -6.84032783e+02 -6.87577029e+02 -6.91121275e+02
  -6.94665521e+02 -6.98209767e+02 -7.01754013e+02 -7.05298259e+02
  -7.08842505e+02 -7.12386751e+02 -7.15930997e+02 -7.19475243e+02
  -7.23019489e+02 -7.26563735e+02 -7.30107981e+02 -7.33652227e+02
  -7.37196473e+02 -7.40740719e+02 -7.44284965e+02 -7.47829211e+02
  -7.51373457e+02 -7.54917703e+02 -7.58461949e+02 -7.62006195e+02
  -7.65550441e+02 -7.69094687e+02 -7.72638933e+02 -7.76183179e+02
  -7.79727424e+02 -7.83271670e+02 -7.86815916e+02 -7.90360162e+02
  -7.93904408e+02 -7.97448654e+02 -8.00992900e+02 -8.04537146e+02
  -8.08081392e+02 -8.11625638e+02 -8.15169884e+02 -8.18714130e+02
  -8.22258376e+02 -8.25802622e+02 -8.29346868e+02 -8.32891114e+02
  -8.36435360e+02 -8.39979606e+02 -8.43523852e+02 -8.47068098e+02
  -8.50612344e+02 -8.54156590e+02 -8.57700836e+02 -8.61245082e+02
  -8.64789328e+02 -8.68333574e+02 -8.71877820e+02 -8.75422066e+02
  -8.78966312e+02 -8.82510558e+02 -8.86054804e+02 -8.89599050e+02
  -8.93143296e+02 -8.96687542e+02 -9.00231788e+02 -9.03776034e+02
  -9.07320280e+02 -9.10864526e+02 -9.14408772e+02 -9.17953018e+02
  -9.21497264e+02 -9.25041509e+02 -9.28585755e+02 -9.32130001e+02
  -9.35674247e+02 -9.39218493e+02 -9.42762739e+02 -9.46306985e+02
  -9.49851231e+02 -9.53395477e+02 -9.56939723e+02 -9.60483969e+02
  -9.64028215e+02 -9.67572461e+02 -9.71116707e+02 -9.74660953e+02
  -9.78205199e+02 -9.81749445e+02 -9.85293691e+02 -9.88837937e+02
  -9.92382183e+02 -9.95926429e+02 -9.99470675e+02 -1.00301492e+03
  -1.00655917e+03 -1.01010341e+03 -1.01364766e+03 -1.01719190e+03
  -1.02073615e+03 -1.02428040e+03 -1.02782464e+03 -1.03136889e+03
  -1.03491313e+03 -1.03845738e+03 -1.04200163e+03 -1.04554587e+03
  -1.04909012e+03 -1.05263436e+03 -1.05617861e+03 -1.05972286e+03
  -1.06326710e+03 -1.06681135e+03 -1.07035559e+03 -1.07389984e+03
  -1.07744409e+03 -1.08098833e+03 -1.08453258e+03 -1.08807682e+03
  -1.09162107e+03 -1.09516532e+03 -1.09870956e+03 -1.10225381e+03
  -1.10579805e+03 -1.10934230e+03 -1.11288655e+03 -1.11643079e+03
  -1.11997504e+03 -1.12351928e+03 -1.12706353e+03 -1.13060778e+03
  -1.13415202e+03 -1.13769627e+03 -1.14124051e+03 -1.14478476e+03
  -1.14832901e+03 -1.15187325e+03 -1.15541750e+03 -1.15896174e+03
  -1.16250599e+03 -1.16605024e+03 -1.16959448e+03 -1.17313873e+03
  -1.17668297e+03 -1.18022722e+03 -1.18377147e+03 -1.18731571e+03
  -1.19085996e+03 -1.19440420e+03 -1.19794845e+03 -1.20149270e+03
  -1.20503694e+03 -1.20858119e+03 -1.21212543e+03 -1.21566968e+03
  -1.21921393e+03 -1.22275817e+03 -1.22630242e+03 -1.22984666e+03
  -1.23339091e+03 -1.23693516e+03 -1.24047940e+03 -1.24402365e+03
  -1.24756789e+03 -1.25111214e+03 -1.25465639e+03 -1.25820063e+03
  -1.26174488e+03 -1.26528912e+03 -1.26883337e+03 -1.27237762e+03
  -1.27592186e+03 -1.27946611e+03 -1.28301035e+03 -1.28655460e+03
  -1.29009884e+03 -1.29364309e+03 -1.29718734e+03 -1.30073158e+03
  -1.30427583e+03 -1.30782007e+03 -1.31136432e+03 -1.31490857e+03
  -1.31845281e+03 -1.32199706e+03 -1.32554130e+03 -1.32908555e+03
  -1.33262980e+03 -1.33617404e+03 -1.33971829e+03 -1.34326253e+03
  -1.34680678e+03 -1.35035103e+03 -1.35389527e+03 -1.35743952e+03
  -1.36098376e+03 -1.36452801e+03 -1.36807226e+03 -1.37161650e+03]]
[-2.82508836 -1.30841001 -0.04797249]
[]
[-5.74678059 -5.41846318 -0.02072321]
[]
[6.16148845 2.89673437 0.10342573]
[]
[ 4.40003501 -3.49550771  0.22898547]
[[6.55085002e-02 6.35935102e+00 1.26531935e+01 1.89470361e+01
  2.52408786e+01 3.15347211e+01 3.78285636e+01 4.41224061e+01
  5.04162487e+01 5.67100912e+01 6.30039337e+01 6.92977762e+01
  7.55916188e+01 8.18854613e+01 8.81793038e+01 9.44731463e+01
  1.00766989e+02 1.07060831e+02 1.13354674e+02 1.19648516e+02
  1.25942359e+02 1.32236201e+02 1.38530044e+02 1.44823886e+02
  1.51117729e+02 1.57411572e+02 1.63705414e+02 1.69999257e+02
  1.76293099e+02 1.82586942e+02 1.88880784e+02 1.95174627e+02
  2.01468469e+02 2.07762312e+02 2.14056154e+02 2.20349997e+02
  2.26643839e+02 2.32937682e+02 2.39231524e+02 2.45525367e+02
  2.51819209e+02 2.58113052e+02 2.64406894e+02 2.70700737e+02
  2.76994579e+02 2.83288422e+02 2.89582264e+02 2.95876107e+02
  3.02169950e+02 3.08463792e+02 3.14757635e+02 3.21051477e+02
  3.27345320e+02 3.33639162e+02 3.39933005e+02 3.46226847e+02
  3.52520690e+02 3.58814532e+02 3.65108375e+02 3.71402217e+02
  3.77696060e+02 3.83989902e+02 3.90283745e+02 3.96577587e+02
  4.02871430e+02 4.09165272e+02 4.15459115e+02 4.21752957e+02
  4.28046800e+02 4.34340642e+02 4.40634485e+02 4.46928328e+02
  4.53222170e+02 4.59516013e+02 4.65809855e+02 4.72103698e+02
  4.78397540e+02 4.84691383e+02 4.90985225e+02 4.97279068e+02
  5.03572910e+02 5.09866753e+02 5.16160595e+02 5.22454438e+02
  5.28748280e+02 5.35042123e+02 5.41335965e+02 5.47629808e+02
  5.53923650e+02 5.60217493e+02 5.66511335e+02 5.72805178e+02
  5.79099020e+02 5.85392863e+02 5.91686705e+02 5.97980548e+02
  6.04274391e+02 6.10568233e+02 6.16862076e+02 6.23155918e+02
  6.29449761e+02 6.35743603e+02 6.42037446e+02 6.48331288e+02
  6.54625131e+02 6.60918973e+02 6.67212816e+02 6.73506658e+02
  6.79800501e+02 6.86094343e+02 6.92388186e+02 6.98682028e+02
  7.04975871e+02 7.11269713e+02 7.17563556e+02 7.23857398e+02
  7.30151241e+02 7.36445083e+02 7.42738926e+02 7.49032769e+02
  7.55326611e+02 7.61620454e+02 7.67914296e+02 7.74208139e+02
  7.80501981e+02 7.86795824e+02 7.93089666e+02 7.99383509e+02
  8.05677351e+02 8.11971194e+02 8.18265036e+02 8.24558879e+02
  8.30852721e+02 8.37146564e+02 8.43440406e+02 8.49734249e+02
  8.56028091e+02 8.62321934e+02 8.68615776e+02 8.74909619e+02
  8.81203461e+02 8.87497304e+02 8.93791147e+02 9.00084989e+02
  9.06378832e+02 9.12672674e+02 9.18966517e+02 9.25260359e+02
  9.31554202e+02 9.37848044e+02 9.44141887e+02 9.50435729e+02
  9.56729572e+02 9.63023414e+02 9.69317257e+02 9.75611099e+02
  9.81904942e+02 9.88198784e+02 9.94492627e+02 1.00078647e+03
  1.00708031e+03 1.01337415e+03 1.01966800e+03 1.02596184e+03
  1.03225568e+03 1.03854952e+03 1.04484337e+03 1.05113721e+03
  1.05743105e+03 1.06372489e+03 1.07001874e+03 1.07631258e+03
  1.08260642e+03 1.08890026e+03 1.09519411e+03 1.10148795e+03
  1.10778179e+03 1.11407563e+03 1.12036948e+03 1.12666332e+03
  1.13295716e+03 1.13925100e+03 1.14554485e+03 1.15183869e+03
  1.15813253e+03 1.16442637e+03 1.17072022e+03 1.17701406e+03
  1.18330790e+03 1.18960175e+03 1.19589559e+03 1.20218943e+03
  1.20848327e+03 1.21477712e+03 1.22107096e+03 1.22736480e+03
  1.23365864e+03 1.23995249e+03 1.24624633e+03 1.25254017e+03
  1.25883401e+03 1.26512786e+03 1.27142170e+03 1.27771554e+03
  1.28400938e+03 1.29030323e+03 1.29659707e+03 1.30289091e+03
  1.30918475e+03 1.31547860e+03 1.32177244e+03 1.32806628e+03
  1.33436012e+03 1.34065397e+03 1.34694781e+03 1.35324165e+03
  1.35953549e+03 1.36582934e+03 1.37212318e+03 1.37841702e+03
  1.38471086e+03 1.39100471e+03 1.39729855e+03 1.40359239e+03
  1.40988623e+03 1.41618008e+03 1.42247392e+03 1.42876776e+03
  1.43506160e+03 1.44135545e+03 1.44764929e+03 1.45394313e+03
  1.46023697e+03 1.46653082e+03 1.47282466e+03 1.47911850e+03
  1.48541234e+03 1.49170619e+03 1.49800003e+03 1.50429387e+03
  1.51058771e+03 1.51688156e+03 1.52317540e+03 1.52946924e+03
  1.53576308e+03 1.54205693e+03 1.54835077e+03 1.55464461e+03
  1.56093845e+03 1.56723230e+03 1.57352614e+03 1.57981998e+03
  1.58611382e+03 1.59240767e+03 1.59870151e+03 1.60499535e+03
  1.61128919e+03 1.61758304e+03 1.62387688e+03 1.63017072e+03
  1.63646456e+03 1.64275841e+03 1.64905225e+03 1.65534609e+03
  1.66163993e+03 1.66793378e+03 1.67422762e+03 1.68052146e+03
  1.68681530e+03 1.69310915e+03 1.69940299e+03 1.70569683e+03
  1.71199067e+03 1.71828452e+03 1.72457836e+03 1.73087220e+03
  1.73716604e+03 1.74345989e+03 1.74975373e+03 1.75604757e+03
  1.76234141e+03 1.76863526e+03 1.77492910e+03 1.78122294e+03
  1.78751678e+03 1.79381063e+03 1.80010447e+03 1.80639831e+03
  1.81269215e+03 1.81898600e+03 1.82527984e+03 1.83157368e+03
  1.83786752e+03 1.84416137e+03 1.85045521e+03 1.85674905e+03
  1.86304289e+03 1.86933674e+03 1.87563058e+03 1.88192442e+03
  1.88821826e+03 1.89451211e+03 1.90080595e+03 1.90709979e+03
  1.91339363e+03 1.91968748e+03 1.92598132e+03 1.93227516e+03
  1.93856901e+03 1.94486285e+03 1.95115669e+03 1.95745053e+03
  1.96374438e+03 1.97003822e+03 1.97633206e+03 1.98262590e+03
  1.98891975e+03 1.99521359e+03 2.00150743e+03 2.00780127e+03
  2.01409512e+03 2.02038896e+03 2.02668280e+03 2.03297664e+03
  2.03927049e+03 2.04556433e+03 2.05185817e+03 2.05815201e+03
  2.06444586e+03 2.07073970e+03 2.07703354e+03 2.08332738e+03
  2.08962123e+03 2.09591507e+03 2.10220891e+03 2.10850275e+03
  2.11479660e+03 2.12109044e+03 2.12738428e+03 2.13367812e+03
  2.13997197e+03 2.14626581e+03 2.15255965e+03 2.15885349e+03
  2.16514734e+03 2.17144118e+03 2.17773502e+03 2.18402886e+03
  2.19032271e+03 2.19661655e+03 2.20291039e+03 2.20920423e+03
  2.21549808e+03 2.22179192e+03 2.22808576e+03 2.23437960e+03
  2.24067345e+03 2.24696729e+03 2.25326113e+03 2.25955497e+03
  2.26584882e+03 2.27214266e+03 2.27843650e+03 2.28473034e+03
  2.29102419e+03 2.29731803e+03 2.30361187e+03 2.30990571e+03
  2.31619956e+03 2.32249340e+03 2.32878724e+03 2.33508108e+03
  2.34137493e+03 2.34766877e+03 2.35396261e+03 2.36025645e+03
  2.36655030e+03 2.37284414e+03 2.37913798e+03 2.38543182e+03
  2.39172567e+03 2.39801951e+03 2.40431335e+03 2.41060719e+03
  2.41690104e+03 2.42319488e+03 2.42948872e+03 2.43578256e+03]]
[-4.19627966 -6.73917435  0.0624475 ]
[[ 9.26634281e-03 -3.10408215e+00 -6.21743064e+00 -9.33077913e+00
  -1.24441276e+01 -1.55574761e+01 -1.86708246e+01 -2.17841731e+01
  -2.48975216e+01 -2.80108701e+01 -3.11242186e+01 -3.42375671e+01
  -3.73509155e+01 -4.04642640e+01 -4.35776125e+01 -4.66909610e+01
  -4.98043095e+01 -5.29176580e+01 -5.60310065e+01 -5.91443550e+01
  -6.22577035e+01 -6.53710520e+01 -6.84844005e+01 -7.15977489e+01
  -7.47110974e+01 -7.78244459e+01 -8.09377944e+01 -8.40511429e+01
  -8.71644914e+01 -9.02778399e+01 -9.33911884e+01 -9.65045369e+01
  -9.96178854e+01 -1.02731234e+02 -1.05844582e+02 -1.08957931e+02
  -1.12071279e+02 -1.15184628e+02 -1.18297976e+02 -1.21411325e+02
  -1.24524673e+02 -1.27638022e+02 -1.30751370e+02 -1.33864719e+02
  -1.36978067e+02 -1.40091416e+02 -1.43204764e+02 -1.46318113e+02
  -1.49431461e+02 -1.52544810e+02 -1.55658158e+02 -1.58771507e+02
  -1.61884855e+02 -1.64998204e+02 -1.68111552e+02 -1.71224901e+02
  -1.74338249e+02 -1.77451598e+02 -1.80564946e+02 -1.83678295e+02
  -1.86791643e+02 -1.89904992e+02 -1.93018340e+02 -1.96131689e+02
  -1.99245037e+02 -2.02358386e+02 -2.05471734e+02 -2.08585083e+02
  -2.11698431e+02 -2.14811780e+02 -2.17925128e+02 -2.21038476e+02
  -2.24151825e+02 -2.27265173e+02 -2.30378522e+02 -2.33491870e+02
  -2.36605219e+02 -2.39718567e+02 -2.42831916e+02 -2.45945264e+02
  -2.49058613e+02 -2.52171961e+02 -2.55285310e+02 -2.58398658e+02
  -2.61512007e+02 -2.64625355e+02 -2.67738704e+02 -2.70852052e+02
  -2.73965401e+02 -2.77078749e+02 -2.80192098e+02 -2.83305446e+02
  -2.86418795e+02 -2.89532143e+02 -2.92645492e+02 -2.95758840e+02
  -2.98872189e+02 -3.01985537e+02 -3.05098886e+02 -3.08212234e+02
  -3.11325583e+02 -3.14438931e+02 -3.17552280e+02 -3.20665628e+02
  -3.23778977e+02 -3.26892325e+02 -3.30005674e+02 -3.33119022e+02
  -3.36232371e+02 -3.39345719e+02 -3.42459068e+02 -3.45572416e+02
  -3.48685765e+02 -3.51799113e+02 -3.54912462e+02 -3.58025810e+02
  -3.61139159e+02 -3.64252507e+02 -3.67365856e+02 -3.70479204e+02
  -3.73592553e+02 -3.76705901e+02 -3.79819250e+02 -3.82932598e+02
  -3.86045946e+02 -3.89159295e+02 -3.92272643e+02 -3.95385992e+02
  -3.98499340e+02 -4.01612689e+02 -4.04726037e+02 -4.07839386e+02
  -4.10952734e+02 -4.14066083e+02 -4.17179431e+02 -4.20292780e+02
  -4.23406128e+02 -4.26519477e+02 -4.29632825e+02 -4.32746174e+02
  -4.35859522e+02 -4.38972871e+02 -4.42086219e+02 -4.45199568e+02
  -4.48312916e+02 -4.51426265e+02 -4.54539613e+02 -4.57652962e+02
  -4.60766310e+02 -4.63879659e+02 -4.66993007e+02 -4.70106356e+02
  -4.73219704e+02 -4.76333053e+02 -4.79446401e+02 -4.82559750e+02
  -4.85673098e+02 -4.88786447e+02 -4.91899795e+02 -4.95013144e+02
  -4.98126492e+02 -5.01239841e+02 -5.04353189e+02 -5.07466538e+02
  -5.10579886e+02 -5.13693235e+02 -5.16806583e+02 -5.19919932e+02
  -5.23033280e+02 -5.26146629e+02 -5.29259977e+02 -5.32373326e+02
  -5.35486674e+02 -5.38600023e+02 -5.41713371e+02 -5.44826720e+02
  -5.47940068e+02 -5.51053416e+02 -5.54166765e+02 -5.57280113e+02
  -5.60393462e+02 -5.63506810e+02 -5.66620159e+02 -5.69733507e+02
  -5.72846856e+02 -5.75960204e+02 -5.79073553e+02 -5.82186901e+02
  -5.85300250e+02 -5.88413598e+02 -5.91526947e+02 -5.94640295e+02
  -5.97753644e+02 -6.00866992e+02 -6.03980341e+02 -6.07093689e+02
  -6.10207038e+02 -6.13320386e+02 -6.16433735e+02 -6.19547083e+02
  -6.22660432e+02 -6.25773780e+02 -6.28887129e+02 -6.32000477e+02
  -6.35113826e+02 -6.38227174e+02 -6.41340523e+02 -6.44453871e+02
  -6.47567220e+02 -6.50680568e+02 -6.53793917e+02 -6.56907265e+02
  -6.60020614e+02 -6.63133962e+02 -6.66247311e+02 -6.69360659e+02
  -6.72474008e+02 -6.75587356e+02 -6.78700705e+02 -6.81814053e+02
  -6.84927402e+02 -6.88040750e+02 -6.91154099e+02 -6.94267447e+02
  -6.97380796e+02 -7.00494144e+02 -7.03607493e+02 -7.06720841e+02
  -7.09834190e+02 -7.12947538e+02 -7.16060887e+02 -7.19174235e+02
  -7.22287583e+02 -7.25400932e+02 -7.28514280e+02 -7.31627629e+02
  -7.34740977e+02 -7.37854326e+02 -7.40967674e+02 -7.44081023e+02
  -7.47194371e+02 -7.50307720e+02 -7.53421068e+02 -7.56534417e+02
  -7.59647765e+02 -7.62761114e+02 -7.65874462e+02 -7.68987811e+02
  -7.72101159e+02 -7.75214508e+02 -7.78327856e+02 -7.81441205e+02
  -7.84554553e+02 -7.87667902e+02 -7.90781250e+02 -7.93894599e+02
  -7.97007947e+02 -8.00121296e+02 -8.03234644e+02 -8.06347993e+02
  -8.09461341e+02 -8.12574690e+02 -8.15688038e+02 -8.18801387e+02
  -8.21914735e+02 -8.25028084e+02 -8.28141432e+02 -8.31254781e+02
  -8.34368129e+02 -8.37481478e+02 -8.40594826e+02 -8.43708175e+02
  -8.46821523e+02 -8.49934872e+02 -8.53048220e+02 -8.56161569e+02
  -8.59274917e+02 -8.62388266e+02 -8.65501614e+02 -8.68614963e+02
  -8.71728311e+02 -8.74841660e+02 -8.77955008e+02 -8.81068357e+02
  -8.84181705e+02 -8.87295053e+02 -8.90408402e+02 -8.93521750e+02
  -8.96635099e+02 -8.99748447e+02 -9.02861796e+02 -9.05975144e+02
  -9.09088493e+02 -9.12201841e+02 -9.15315190e+02 -9.18428538e+02
  -9.21541887e+02 -9.24655235e+02 -9.27768584e+02 -9.30881932e+02
  -9.33995281e+02 -9.37108629e+02 -9.40221978e+02 -9.43335326e+02
  -9.46448675e+02 -9.49562023e+02 -9.52675372e+02 -9.55788720e+02
  -9.58902069e+02 -9.62015417e+02 -9.65128766e+02 -9.68242114e+02
  -9.71355463e+02 -9.74468811e+02 -9.77582160e+02 -9.80695508e+02
  -9.83808857e+02 -9.86922205e+02 -9.90035554e+02 -9.93148902e+02
  -9.96262251e+02 -9.99375599e+02 -1.00248895e+03 -1.00560230e+03
  -1.00871564e+03 -1.01182899e+03 -1.01494234e+03 -1.01805569e+03
  -1.02116904e+03 -1.02428239e+03 -1.02739574e+03 -1.03050908e+03
  -1.03362243e+03 -1.03673578e+03 -1.03984913e+03 -1.04296248e+03
  -1.04607583e+03 -1.04918917e+03 -1.05230252e+03 -1.05541587e+03
  -1.05852922e+03 -1.06164257e+03 -1.06475592e+03 -1.06786927e+03
  -1.07098261e+03 -1.07409596e+03 -1.07720931e+03 -1.08032266e+03
  -1.08343601e+03 -1.08654936e+03 -1.08966271e+03 -1.09277605e+03
  -1.09588940e+03 -1.09900275e+03 -1.10211610e+03 -1.10522945e+03
  -1.10834280e+03 -1.11145614e+03 -1.11456949e+03 -1.11768284e+03
  -1.12079619e+03 -1.12390954e+03 -1.12702289e+03 -1.13013624e+03
  -1.13324958e+03 -1.13636293e+03 -1.13947628e+03 -1.14258963e+03
  -1.14570298e+03 -1.14881633e+03 -1.15192968e+03 -1.15504302e+03
  -1.15815637e+03 -1.16126972e+03 -1.16438307e+03 -1.16749642e+03
  -1.17060977e+03 -1.17372311e+03 -1.17683646e+03 -1.17994981e+03
  -1.18306316e+03 -1.18617651e+03 -1.18928986e+03 -1.19240321e+03
  -1.19551655e+03 -1.19862990e+03 -1.20174325e+03 -1.20485660e+03]]
[-4.22792254 -5.40789749  0.02438571]
[[ 4.50927688e-03 -3.90451687e+00 -7.81354301e+00 -1.17225692e+01
  -1.56315953e+01 -1.95406214e+01 -2.34496476e+01 -2.73586737e+01
  -3.12676999e+01 -3.51767260e+01 -3.90857522e+01 -4.29947783e+01
  -4.69038045e+01 -5.08128306e+01 -5.47218568e+01 -5.86308829e+01
  -6.25399090e+01 -6.64489352e+01 -7.03579613e+01 -7.42669875e+01
  -7.81760136e+01 -8.20850398e+01 -8.59940659e+01 -8.99030921e+01
  -9.38121182e+01 -9.77211444e+01 -1.01630170e+02 -1.05539197e+02
  -1.09448223e+02 -1.13357249e+02 -1.17266275e+02 -1.21175301e+02
  -1.25084327e+02 -1.28993354e+02 -1.32902380e+02 -1.36811406e+02
  -1.40720432e+02 -1.44629458e+02 -1.48538484e+02 -1.52447510e+02
  -1.56356537e+02 -1.60265563e+02 -1.64174589e+02 -1.68083615e+02
  -1.71992641e+02 -1.75901667e+02 -1.79810693e+02 -1.83719720e+02
  -1.87628746e+02 -1.91537772e+02 -1.95446798e+02 -1.99355824e+02
  -2.03264850e+02 -2.07173876e+02 -2.11082903e+02 -2.14991929e+02
  -2.18900955e+02 -2.22809981e+02 -2.26719007e+02 -2.30628033e+02
  -2.34537059e+02 -2.38446086e+02 -2.42355112e+02 -2.46264138e+02
  -2.50173164e+02 -2.54082190e+02 -2.57991216e+02 -2.61900242e+02
  -2.65809269e+02 -2.69718295e+02 -2.73627321e+02 -2.77536347e+02
  -2.81445373e+02 -2.85354399e+02 -2.89263425e+02 -2.93172452e+02
  -2.97081478e+02 -3.00990504e+02 -3.04899530e+02 -3.08808556e+02
  -3.12717582e+02 -3.16626608e+02 -3.20535635e+02 -3.24444661e+02
  -3.28353687e+02 -3.32262713e+02 -3.36171739e+02 -3.40080765e+02
  -3.43989792e+02 -3.47898818e+02 -3.51807844e+02 -3.55716870e+02
  -3.59625896e+02 -3.63534922e+02 -3.67443948e+02 -3.71352975e+02
  -3.75262001e+02 -3.79171027e+02 -3.83080053e+02 -3.86989079e+02
  -3.90898105e+02 -3.94807131e+02 -3.98716158e+02 -4.02625184e+02
  -4.06534210e+02 -4.10443236e+02 -4.14352262e+02 -4.18261288e+02
  -4.22170314e+02 -4.26079341e+02 -4.29988367e+02 -4.33897393e+02
  -4.37806419e+02 -4.41715445e+02 -4.45624471e+02 -4.49533497e+02
  -4.53442524e+02 -4.57351550e+02 -4.61260576e+02 -4.65169602e+02
  -4.69078628e+02 -4.72987654e+02 -4.76896680e+02 -4.80805707e+02
  -4.84714733e+02 -4.88623759e+02 -4.92532785e+02 -4.96441811e+02
  -5.00350837e+02 -5.04259863e+02 -5.08168890e+02 -5.12077916e+02
  -5.15986942e+02 -5.19895968e+02 -5.23804994e+02 -5.27714020e+02
  -5.31623046e+02 -5.35532073e+02 -5.39441099e+02 -5.43350125e+02
  -5.47259151e+02 -5.51168177e+02 -5.55077203e+02 -5.58986229e+02
  -5.62895256e+02 -5.66804282e+02 -5.70713308e+02 -5.74622334e+02
  -5.78531360e+02 -5.82440386e+02 -5.86349413e+02 -5.90258439e+02
  -5.94167465e+02 -5.98076491e+02 -6.01985517e+02 -6.05894543e+02
  -6.09803569e+02 -6.13712596e+02 -6.17621622e+02 -6.21530648e+02
  -6.25439674e+02 -6.29348700e+02 -6.33257726e+02 -6.37166752e+02
  -6.41075779e+02 -6.44984805e+02 -6.48893831e+02 -6.52802857e+02
  -6.56711883e+02 -6.60620909e+02 -6.64529935e+02 -6.68438962e+02
  -6.72347988e+02 -6.76257014e+02 -6.80166040e+02 -6.84075066e+02
  -6.87984092e+02 -6.91893118e+02 -6.95802145e+02 -6.99711171e+02
  -7.03620197e+02 -7.07529223e+02 -7.11438249e+02 -7.15347275e+02
  -7.19256301e+02 -7.23165328e+02 -7.27074354e+02 -7.30983380e+02
  -7.34892406e+02 -7.38801432e+02 -7.42710458e+02 -7.46619484e+02
  -7.50528511e+02 -7.54437537e+02 -7.58346563e+02 -7.62255589e+02
  -7.66164615e+02 -7.70073641e+02 -7.73982667e+02 -7.77891694e+02
  -7.81800720e+02 -7.85709746e+02 -7.89618772e+02 -7.93527798e+02
  -7.97436824e+02 -8.01345850e+02 -8.05254877e+02 -8.09163903e+02
  -8.13072929e+02 -8.16981955e+02 -8.20890981e+02 -8.24800007e+02
  -8.28709034e+02 -8.32618060e+02 -8.36527086e+02 -8.40436112e+02
  -8.44345138e+02 -8.48254164e+02 -8.52163190e+02 -8.56072217e+02
  -8.59981243e+02 -8.63890269e+02 -8.67799295e+02 -8.71708321e+02
  -8.75617347e+02 -8.79526373e+02 -8.83435400e+02 -8.87344426e+02
  -8.91253452e+02 -8.95162478e+02 -8.99071504e+02 -9.02980530e+02
  -9.06889556e+02 -9.10798583e+02 -9.14707609e+02 -9.18616635e+02
  -9.22525661e+02 -9.26434687e+02 -9.30343713e+02 -9.34252739e+02
  -9.38161766e+02 -9.42070792e+02 -9.45979818e+02 -9.49888844e+02
  -9.53797870e+02 -9.57706896e+02 -9.61615922e+02 -9.65524949e+02
  -9.69433975e+02 -9.73343001e+02 -9.77252027e+02 -9.81161053e+02
  -9.85070079e+02 -9.88979105e+02 -9.92888132e+02 -9.96797158e+02
  -1.00070618e+03 -1.00461521e+03 -1.00852424e+03 -1.01243326e+03
  -1.01634229e+03 -1.02025131e+03 -1.02416034e+03 -1.02806937e+03
  -1.03197839e+03 -1.03588742e+03 -1.03979645e+03 -1.04370547e+03
  -1.04761450e+03 -1.05152352e+03 -1.05543255e+03 -1.05934158e+03
  -1.06325060e+03 -1.06715963e+03 -1.07106865e+03 -1.07497768e+03
  -1.07888671e+03 -1.08279573e+03 -1.08670476e+03 -1.09061379e+03
  -1.09452281e+03 -1.09843184e+03 -1.10234086e+03 -1.10624989e+03
  -1.11015892e+03 -1.11406794e+03 -1.11797697e+03 -1.12188599e+03
  -1.12579502e+03 -1.12970405e+03 -1.13361307e+03 -1.13752210e+03
  -1.14143113e+03 -1.14534015e+03 -1.14924918e+03 -1.15315820e+03
  -1.15706723e+03 -1.16097626e+03 -1.16488528e+03 -1.16879431e+03
  -1.17270333e+03 -1.17661236e+03 -1.18052139e+03 -1.18443041e+03
  -1.18833944e+03 -1.19224847e+03 -1.19615749e+03 -1.20006652e+03
  -1.20397554e+03 -1.20788457e+03 -1.21179360e+03 -1.21570262e+03
  -1.21961165e+03 -1.22352067e+03 -1.22742970e+03 -1.23133873e+03
  -1.23524775e+03 -1.23915678e+03 -1.24306580e+03 -1.24697483e+03
  -1.25088386e+03 -1.25479288e+03 -1.25870191e+03 -1.26261094e+03
  -1.26651996e+03 -1.27042899e+03 -1.27433801e+03 -1.27824704e+03
  -1.28215607e+03 -1.28606509e+03 -1.28997412e+03 -1.29388314e+03
  -1.29779217e+03 -1.30170120e+03 -1.30561022e+03 -1.30951925e+03
  -1.31342828e+03 -1.31733730e+03 -1.32124633e+03 -1.32515535e+03
  -1.32906438e+03 -1.33297341e+03 -1.33688243e+03 -1.34079146e+03
  -1.34470048e+03 -1.34860951e+03 -1.35251854e+03 -1.35642756e+03
  -1.36033659e+03 -1.36424562e+03 -1.36815464e+03 -1.37206367e+03
  -1.37597269e+03 -1.37988172e+03 -1.38379075e+03 -1.38769977e+03
  -1.39160880e+03 -1.39551782e+03 -1.39942685e+03 -1.40333588e+03
  -1.40724490e+03 -1.41115393e+03 -1.41506296e+03 -1.41897198e+03
  -1.42288101e+03 -1.42679003e+03 -1.43069906e+03 -1.43460809e+03
  -1.43851711e+03 -1.44242614e+03 -1.44633516e+03 -1.45024419e+03
  -1.45415322e+03 -1.45806224e+03 -1.46197127e+03 -1.46588030e+03
  -1.46978932e+03 -1.47369835e+03 -1.47760737e+03 -1.48151640e+03
  -1.48542543e+03 -1.48933445e+03 -1.49324348e+03 -1.49715250e+03
  -1.50106153e+03 -1.50497056e+03 -1.50887958e+03 -1.51278861e+03]]
[ 3.77772512 -6.25315264  0.28726556]
[[4.59393165e-02 3.06659572e+00 6.08725213e+00 9.10790853e+00
  1.21285649e+01 1.51492213e+01 1.81698778e+01 2.11905342e+01
  2.42111906e+01 2.72318470e+01 3.02525034e+01 3.32731598e+01
  3.62938162e+01 3.93144726e+01 4.23351290e+01 4.53557854e+01
  4.83764418e+01 5.13970982e+01 5.44177546e+01 5.74384110e+01
  6.04590674e+01 6.34797238e+01 6.65003803e+01 6.95210367e+01
  7.25416931e+01 7.55623495e+01 7.85830059e+01 8.16036623e+01
  8.46243187e+01 8.76449751e+01 9.06656315e+01 9.36862879e+01
  9.67069443e+01 9.97276007e+01 1.02748257e+02 1.05768914e+02
  1.08789570e+02 1.11810226e+02 1.14830883e+02 1.17851539e+02
  1.20872196e+02 1.23892852e+02 1.26913508e+02 1.29934165e+02
  1.32954821e+02 1.35975478e+02 1.38996134e+02 1.42016790e+02
  1.45037447e+02 1.48058103e+02 1.51078760e+02 1.54099416e+02
  1.57120072e+02 1.60140729e+02 1.63161385e+02 1.66182042e+02
  1.69202698e+02 1.72223354e+02 1.75244011e+02 1.78264667e+02
  1.81285324e+02 1.84305980e+02 1.87326636e+02 1.90347293e+02
  1.93367949e+02 1.96388606e+02 1.99409262e+02 2.02429919e+02
  2.05450575e+02 2.08471231e+02 2.11491888e+02 2.14512544e+02
  2.17533201e+02 2.20553857e+02 2.23574513e+02 2.26595170e+02
  2.29615826e+02 2.32636483e+02 2.35657139e+02 2.38677795e+02
  2.41698452e+02 2.44719108e+02 2.47739765e+02 2.50760421e+02
  2.53781077e+02 2.56801734e+02 2.59822390e+02 2.62843047e+02
  2.65863703e+02 2.68884359e+02 2.71905016e+02 2.74925672e+02
  2.77946329e+02 2.80966985e+02 2.83987641e+02 2.87008298e+02
  2.90028954e+02 2.93049611e+02 2.96070267e+02 2.99090924e+02
  3.02111580e+02 3.05132236e+02 3.08152893e+02 3.11173549e+02
  3.14194206e+02 3.17214862e+02 3.20235518e+02 3.23256175e+02
  3.26276831e+02 3.29297488e+02 3.32318144e+02 3.35338800e+02
  3.38359457e+02 3.41380113e+02 3.44400770e+02 3.47421426e+02
  3.50442082e+02 3.53462739e+02 3.56483395e+02 3.59504052e+02
  3.62524708e+02 3.65545364e+02 3.68566021e+02 3.71586677e+02
  3.74607334e+02 3.77627990e+02 3.80648646e+02 3.83669303e+02
  3.86689959e+02 3.89710616e+02 3.92731272e+02 3.95751929e+02
  3.98772585e+02 4.01793241e+02 4.04813898e+02 4.07834554e+02
  4.10855211e+02 4.13875867e+02 4.16896523e+02 4.19917180e+02
  4.22937836e+02 4.25958493e+02 4.28979149e+02 4.31999805e+02
  4.35020462e+02 4.38041118e+02 4.41061775e+02 4.44082431e+02
  4.47103087e+02 4.50123744e+02 4.53144400e+02 4.56165057e+02
  4.59185713e+02 4.62206369e+02 4.65227026e+02 4.68247682e+02
  4.71268339e+02 4.74288995e+02 4.77309651e+02 4.80330308e+02
  4.83350964e+02 4.86371621e+02 4.89392277e+02 4.92412934e+02
  4.95433590e+02 4.98454246e+02 5.01474903e+02 5.04495559e+02
  5.07516216e+02 5.10536872e+02 5.13557528e+02 5.16578185e+02
  5.19598841e+02 5.22619498e+02 5.25640154e+02 5.28660810e+02
  5.31681467e+02 5.34702123e+02 5.37722780e+02 5.40743436e+02
  5.43764092e+02 5.46784749e+02 5.49805405e+02 5.52826062e+02
  5.55846718e+02 5.58867374e+02 5.61888031e+02 5.64908687e+02
  5.67929344e+02 5.70950000e+02 5.73970656e+02 5.76991313e+02
  5.80011969e+02 5.83032626e+02 5.86053282e+02 5.89073939e+02
  5.92094595e+02 5.95115251e+02 5.98135908e+02 6.01156564e+02
  6.04177221e+02 6.07197877e+02 6.10218533e+02 6.13239190e+02
  6.16259846e+02 6.19280503e+02 6.22301159e+02 6.25321815e+02
  6.28342472e+02 6.31363128e+02 6.34383785e+02 6.37404441e+02
  6.40425097e+02 6.43445754e+02 6.46466410e+02 6.49487067e+02
  6.52507723e+02 6.55528379e+02 6.58549036e+02 6.61569692e+02
  6.64590349e+02 6.67611005e+02 6.70631661e+02 6.73652318e+02
  6.76672974e+02 6.79693631e+02 6.82714287e+02 6.85734944e+02
  6.88755600e+02 6.91776256e+02 6.94796913e+02 6.97817569e+02
  7.00838226e+02 7.03858882e+02 7.06879538e+02 7.09900195e+02
  7.12920851e+02 7.15941508e+02 7.18962164e+02 7.21982820e+02
  7.25003477e+02 7.28024133e+02 7.31044790e+02 7.34065446e+02
  7.37086102e+02 7.40106759e+02 7.43127415e+02 7.46148072e+02
  7.49168728e+02 7.52189384e+02 7.55210041e+02 7.58230697e+02
  7.61251354e+02 7.64272010e+02 7.67292666e+02 7.70313323e+02
  7.73333979e+02 7.76354636e+02 7.79375292e+02 7.82395949e+02
  7.85416605e+02 7.88437261e+02 7.91457918e+02 7.94478574e+02
  7.97499231e+02 8.00519887e+02 8.03540543e+02 8.06561200e+02
  8.09581856e+02 8.12602513e+02 8.15623169e+02 8.18643825e+02
  8.21664482e+02 8.24685138e+02 8.27705795e+02 8.30726451e+02
  8.33747107e+02 8.36767764e+02 8.39788420e+02 8.42809077e+02
  8.45829733e+02 8.48850389e+02 8.51871046e+02 8.54891702e+02
  8.57912359e+02 8.60933015e+02 8.63953671e+02 8.66974328e+02
  8.69994984e+02 8.73015641e+02 8.76036297e+02 8.79056954e+02
  8.82077610e+02 8.85098266e+02 8.88118923e+02 8.91139579e+02
  8.94160236e+02 8.97180892e+02 9.00201548e+02 9.03222205e+02
  9.06242861e+02 9.09263518e+02 9.12284174e+02 9.15304830e+02
  9.18325487e+02 9.21346143e+02 9.24366800e+02 9.27387456e+02
  9.30408112e+02 9.33428769e+02 9.36449425e+02 9.39470082e+02
  9.42490738e+02 9.45511394e+02 9.48532051e+02 9.51552707e+02
  9.54573364e+02 9.57594020e+02 9.60614676e+02 9.63635333e+02
  9.66655989e+02 9.69676646e+02 9.72697302e+02 9.75717958e+02
  9.78738615e+02 9.81759271e+02 9.84779928e+02 9.87800584e+02
  9.90821241e+02 9.93841897e+02 9.96862553e+02 9.99883210e+02
  1.00290387e+03 1.00592452e+03 1.00894518e+03 1.01196584e+03
  1.01498649e+03 1.01800715e+03 1.02102780e+03 1.02404846e+03
  1.02706912e+03 1.03008977e+03 1.03311043e+03 1.03613109e+03
  1.03915174e+03 1.04217240e+03 1.04519306e+03 1.04821371e+03
  1.05123437e+03 1.05425503e+03 1.05727568e+03 1.06029634e+03
  1.06331699e+03 1.06633765e+03 1.06935831e+03 1.07237896e+03
  1.07539962e+03 1.07842028e+03 1.08144093e+03 1.08446159e+03
  1.08748225e+03 1.09050290e+03 1.09352356e+03 1.09654421e+03
  1.09956487e+03 1.10258553e+03 1.10560618e+03 1.10862684e+03
  1.11164750e+03 1.11466815e+03 1.11768881e+03 1.12070947e+03
  1.12373012e+03 1.12675078e+03 1.12977144e+03 1.13279209e+03
  1.13581275e+03 1.13883340e+03 1.14185406e+03 1.14487472e+03
  1.14789537e+03 1.15091603e+03 1.15393669e+03 1.15695734e+03
  1.15997800e+03 1.16299866e+03 1.16601931e+03 1.16903997e+03]]
[ 0.9826044   5.6543793  -0.12827046]
[[ 2.26851522e-02 -8.46202787e-01 -1.71509073e+00 -2.58397866e+00
  -3.45286660e+00 -4.32175454e+00 -5.19064248e+00 -6.05953042e+00
  -6.92841836e+00 -7.79730630e+00 -8.66619424e+00 -9.53508218e+00
  -1.04039701e+01 -1.12728581e+01 -1.21417460e+01 -1.30106339e+01
  -1.38795219e+01 -1.47484098e+01 -1.56172978e+01 -1.64861857e+01
  -1.73550736e+01 -1.82239616e+01 -1.90928495e+01 -1.99617374e+01
  -2.08306254e+01 -2.16995133e+01 -2.25684013e+01 -2.34372892e+01
  -2.43061771e+01 -2.51750651e+01 -2.60439530e+01 -2.69128410e+01
  -2.77817289e+01 -2.86506168e+01 -2.95195048e+01 -3.03883927e+01
  -3.12572807e+01 -3.21261686e+01 -3.29950565e+01 -3.38639445e+01
  -3.47328324e+01 -3.56017203e+01 -3.64706083e+01 -3.73394962e+01
  -3.82083842e+01 -3.90772721e+01 -3.99461600e+01 -4.08150480e+01
  -4.16839359e+01 -4.25528239e+01 -4.34217118e+01 -4.42905997e+01
  -4.51594877e+01 -4.60283756e+01 -4.68972636e+01 -4.77661515e+01
  -4.86350394e+01 -4.95039274e+01 -5.03728153e+01 -5.12417033e+01
  -5.21105912e+01 -5.29794791e+01 -5.38483671e+01 -5.47172550e+01
  -5.55861429e+01 -5.64550309e+01 -5.73239188e+01 -5.81928068e+01
  -5.90616947e+01 -5.99305826e+01 -6.07994706e+01 -6.16683585e+01
  -6.25372465e+01 -6.34061344e+01 -6.42750223e+01 -6.51439103e+01
  -6.60127982e+01 -6.68816862e+01 -6.77505741e+01 -6.86194620e+01
  -6.94883500e+01 -7.03572379e+01 -7.12261258e+01 -7.20950138e+01
  -7.29639017e+01 -7.38327897e+01 -7.47016776e+01 -7.55705655e+01
  -7.64394535e+01 -7.73083414e+01 -7.81772294e+01 -7.90461173e+01
  -7.99150052e+01 -8.07838932e+01 -8.16527811e+01 -8.25216691e+01
  -8.33905570e+01 -8.42594449e+01 -8.51283329e+01 -8.59972208e+01
  -8.68661088e+01 -8.77349967e+01 -8.86038846e+01 -8.94727726e+01
  -9.03416605e+01 -9.12105484e+01 -9.20794364e+01 -9.29483243e+01
  -9.38172123e+01 -9.46861002e+01 -9.55549881e+01 -9.64238761e+01
  -9.72927640e+01 -9.81616520e+01 -9.90305399e+01 -9.98994278e+01
  -1.00768316e+02 -1.01637204e+02 -1.02506092e+02 -1.03374980e+02
  -1.04243868e+02 -1.05112755e+02 -1.05981643e+02 -1.06850531e+02
  -1.07719419e+02 -1.08588307e+02 -1.09457195e+02 -1.10326083e+02
  -1.11194971e+02 -1.12063859e+02 -1.12932747e+02 -1.13801635e+02
  -1.14670523e+02 -1.15539411e+02 -1.16408299e+02 -1.17277187e+02
  -1.18146075e+02 -1.19014962e+02 -1.19883850e+02 -1.20752738e+02
  -1.21621626e+02 -1.22490514e+02 -1.23359402e+02 -1.24228290e+02
  -1.25097178e+02 -1.25966066e+02 -1.26834954e+02 -1.27703842e+02
  -1.28572730e+02 -1.29441618e+02 -1.30310506e+02 -1.31179394e+02
  -1.32048282e+02 -1.32917170e+02 -1.33786057e+02 -1.34654945e+02
  -1.35523833e+02 -1.36392721e+02 -1.37261609e+02 -1.38130497e+02
  -1.38999385e+02 -1.39868273e+02 -1.40737161e+02 -1.41606049e+02
  -1.42474937e+02 -1.43343825e+02 -1.44212713e+02 -1.45081601e+02
  -1.45950489e+02 -1.46819377e+02 -1.47688264e+02 -1.48557152e+02
  -1.49426040e+02 -1.50294928e+02 -1.51163816e+02 -1.52032704e+02
  -1.52901592e+02 -1.53770480e+02 -1.54639368e+02 -1.55508256e+02
  -1.56377144e+02 -1.57246032e+02 -1.58114920e+02 -1.58983808e+02
  -1.59852696e+02 -1.60721584e+02 -1.61590472e+02 -1.62459359e+02
  -1.63328247e+02 -1.64197135e+02 -1.65066023e+02 -1.65934911e+02
  -1.66803799e+02 -1.67672687e+02 -1.68541575e+02 -1.69410463e+02
  -1.70279351e+02 -1.71148239e+02 -1.72017127e+02 -1.72886015e+02
  -1.73754903e+02 -1.74623791e+02 -1.75492679e+02 -1.76361566e+02
  -1.77230454e+02 -1.78099342e+02 -1.78968230e+02 -1.79837118e+02
  -1.80706006e+02 -1.81574894e+02 -1.82443782e+02 -1.83312670e+02
  -1.84181558e+02 -1.85050446e+02 -1.85919334e+02 -1.86788222e+02
  -1.87657110e+02 -1.88525998e+02 -1.89394886e+02 -1.90263774e+02
  -1.91132661e+02 -1.92001549e+02 -1.92870437e+02 -1.93739325e+02
  -1.94608213e+02 -1.95477101e+02 -1.96345989e+02 -1.97214877e+02
  -1.98083765e+02 -1.98952653e+02 -1.99821541e+02 -2.00690429e+02
  -2.01559317e+02 -2.02428205e+02 -2.03297093e+02 -2.04165981e+02
  -2.05034868e+02 -2.05903756e+02 -2.06772644e+02 -2.07641532e+02
  -2.08510420e+02 -2.09379308e+02 -2.10248196e+02 -2.11117084e+02
  -2.11985972e+02 -2.12854860e+02 -2.13723748e+02 -2.14592636e+02
  -2.15461524e+02 -2.16330412e+02 -2.17199300e+02 -2.18068188e+02
  -2.18937075e+02 -2.19805963e+02 -2.20674851e+02 -2.21543739e+02
  -2.22412627e+02 -2.23281515e+02 -2.24150403e+02 -2.25019291e+02
  -2.25888179e+02 -2.26757067e+02 -2.27625955e+02 -2.28494843e+02
  -2.29363731e+02 -2.30232619e+02 -2.31101507e+02 -2.31970395e+02
  -2.32839283e+02 -2.33708170e+02 -2.34577058e+02 -2.35445946e+02
  -2.36314834e+02 -2.37183722e+02 -2.38052610e+02 -2.38921498e+02
  -2.39790386e+02 -2.40659274e+02 -2.41528162e+02 -2.42397050e+02
  -2.43265938e+02 -2.44134826e+02 -2.45003714e+02 -2.45872602e+02
  -2.46741490e+02 -2.47610377e+02 -2.48479265e+02 -2.49348153e+02
  -2.50217041e+02 -2.51085929e+02 -2.51954817e+02 -2.52823705e+02
  -2.53692593e+02 -2.54561481e+02 -2.55430369e+02 -2.56299257e+02
  -2.57168145e+02 -2.58037033e+02 -2.58905921e+02 -2.59774809e+02
  -2.60643697e+02 -2.61512585e+02 -2.62381472e+02 -2.63250360e+02
  -2.64119248e+02 -2.64988136e+02 -2.65857024e+02 -2.66725912e+02
  -2.67594800e+02 -2.68463688e+02 -2.69332576e+02 -2.70201464e+02
  -2.71070352e+02 -2.71939240e+02 -2.72808128e+02 -2.73677016e+02
  -2.74545904e+02 -2.75414792e+02 -2.76283679e+02 -2.77152567e+02
  -2.78021455e+02 -2.78890343e+02 -2.79759231e+02 -2.80628119e+02
  -2.81497007e+02 -2.82365895e+02 -2.83234783e+02 -2.84103671e+02
  -2.84972559e+02 -2.85841447e+02 -2.86710335e+02 -2.87579223e+02
  -2.88448111e+02 -2.89316999e+02 -2.90185886e+02 -2.91054774e+02
  -2.91923662e+02 -2.92792550e+02 -2.93661438e+02 -2.94530326e+02
  -2.95399214e+02 -2.96268102e+02 -2.97136990e+02 -2.98005878e+02
  -2.98874766e+02 -2.99743654e+02 -3.00612542e+02 -3.01481430e+02
  -3.02350318e+02 -3.03219206e+02 -3.04088094e+02 -3.04956981e+02
  -3.05825869e+02 -3.06694757e+02 -3.07563645e+02 -3.08432533e+02
  -3.09301421e+02 -3.10170309e+02 -3.11039197e+02 -3.11908085e+02
  -3.12776973e+02 -3.13645861e+02 -3.14514749e+02 -3.15383637e+02
  -3.16252525e+02 -3.17121413e+02 -3.17990301e+02 -3.18859188e+02
  -3.19728076e+02 -3.20596964e+02 -3.21465852e+02 -3.22334740e+02
  -3.23203628e+02 -3.24072516e+02 -3.24941404e+02 -3.25810292e+02
  -3.26679180e+02 -3.27548068e+02 -3.28416956e+02 -3.29285844e+02
  -3.30154732e+02 -3.31023620e+02 -3.31892508e+02 -3.32761396e+02
  -3.33630283e+02 -3.34499171e+02 -3.35368059e+02 -3.36236947e+02]]

Step 7: Camera rotation and translation (four solutions)

Factorize essential matrix $E=[T]_x R$ where $R$ is rotation and $T$ is a translation. Find solutions $R_1$, $R_2$ and $T_1$, $T_2$. Use camera 1 for world coordinates. Define projection matrix for camera 1 as $P_w = [I|0]$ and compute four projection matrices for the second camera $P_a$, $P_b$, $P_c$, $P_d$.

Hint 1: for array multiplication use $dot$ or $matmul$, never $*$.

Hint 2: function $svd$ from $linalg$ returns $V^T$ rather than $V$ (the 2nd orthogonal matrix in svd decomposition $E = USV^T$).

Warning: remember that python uses 0 as a starting index for the rows or columns in arrays. For example, $A[0]$ denotes the first row of matrix $A$, while $P_w[2]$ stands for the 3rd row of the corresponding projection matrix and $E[:,[1]]$ is the second column of the essential matrix.

W = np.reshape(a=[0,-1, 0, 1, 0, 0, 0, 0, 1], newshape=(3,3))
R1 = np.matmul(np.matmul(Ue, W), Ve)
R2 = np.matmul(np.matmul(Ue, W.T), Ve)
T1 = Ue[:, -1]
T2 = -Ue[:, -1]

# first camera matrix
Pw = np.column_stack((np.identity(3), np.zeros(3)))

# four possible matrices for the second camera
Pa = np.column_stack((R1, T1))
Pb = np.column_stack((R1, T2))
Pc = np.column_stack((R2, T1))
Pd = np.column_stack((R2, T2))

Summary of Structure-from-Motion (the remaining steps 8-11):

In these 3D reconstruction steps you should use the world coordinate system consistent with the projection matrices estimated in step 7. In all steps you should obtain solutions for all four distinct cases of the second camera: $P_a$, $P_b$, $P_c$, $P_d$. First, step 8 is to Implement least squares (you can use $svd$ or $inv$ functions) for "triangulating" 3D points corresponding to pairs of matched features that are inliers for estimated $E$ (i.e. consistent with the epipolar geometry). Make sure to use normalized coordinates for image points. Then (step 9) you will compute camera positioning (optical centers and calibrated image centers as 3D points) in the world coordinate system. This is used in data visualsation step 10 (fully implemented). That step visulaizes in 3D both camera positions (red - optcal centers, green - image centers) and triangulated points (blue) for four possible cases of the second camera. You should identify one case when solution has 3D points in front of both cameras. In the last step 11 you will project 3D points onto each camera, convert to uncalibrated coordinates, and display these projected points (use red) together with the original features (use blue). Observe if the are red and blue points are close in each image.

Step 8: Triangulation (four solutions)

# Select normalized coordinates for matched features that are inliers for essential matrix E. 
# Form matrix A in equation AX=0 where X represent 4 vectors (homogeneous representation of 3D point).
# Use your solution for Problem 6.
# Each camera (projection matrix P) will define its own A  

# HINT: to keep it simple, first solve the problem for one match.

n_ptL = n_ptsL[ind_sample, :][0]
n_ptR = n_ptsR[ind_sample, :][0]

def construct_A(pA, pB, xA, xB):
    A = np.reshape(a=[xA[1]*pA[2, :] -pA[1, :], pA[0, :] -xA[0]*pA[2, :], 
                      xB[1]*pB[2, :] -pB[1, :], pB[0, :] -xB[0]*pB[2, :]],
                      newshape=(4, 4))
    return A

# print(construct_A(Pw, Pa, n_ptL, n_ptR))

Aa = construct_A(Pw, Pa, n_ptL, n_ptR)

Ab = construct_A(Pw, Pb, n_ptL, n_ptR)

Ac = construct_A(Pw, Pc, n_ptL, n_ptR)

Ad = construct_A(Pw, Pd, n_ptL, n_ptR)

Solution using least squares: assume homogeneous 3D point $X=[X_1,X_2,X_3,1]$. Then, $AX=0$ gives 4 equations for 3 unknowns. Use approach 1 (inhomogeneous least squares) discussed for homography estimation (Topic 6).

Xa = []
Xb = []
Xc = []
Xd = []

n_ptL = n_ptsL
n_ptR = n_ptsR
for i, j in EmatchesRansac:

    Aa = construct_A(Pw, Pa, n_ptL[i], n_ptR[j])
    Ab = construct_A(Pw, Pb, n_ptL[i], n_ptR[j])
    Ac = construct_A(Pw, Pc, n_ptL[i], n_ptR[j])
    Ad = construct_A(Pw, Pd, n_ptL[i], n_ptR[j])

    # least squares for solving linear system A_{0:2} X_{0:2} = - A_3 
    Aa_02 = Aa[:, :3]       # the first 3 columns of 4x4 matrix A
    Aa_3  = Aa[:, 3].reshape(4, 1)       # the last column on 4x4 matrix A
    Ab_02 = Ab[:, :3]    
    Ab_3  = Ab[:, 3]
    Ac_02 = Ac[:, :3]
    Ac_3  = Ac[:, 3]
    Ad_02 = Ad[:, :3]
    Ad_3  = Ad[:, 3]


    # Nx3 matrices: N rows with 3D point coordinates for N reconstructed points (N=num_inliers)

    Xa.append(np.matmul(np.matmul(np.linalg.inv(np.matmul(Aa_02.T, Aa_02)), Aa_02.T), Aa_3))
    Xb.append(np.matmul(np.matmul(np.linalg.inv(np.matmul(Ab_02.T, Ab_02)), Ab_02.T), Ab_3))
    Xc.append(np.matmul(np.matmul(np.linalg.inv(np.matmul(Ac_02.T, Ac_02)), Ac_02.T), Ac_3))
    Xd.append(np.matmul(np.matmul(np.linalg.inv(np.matmul(Ad_02.T, Ad_02)), Ad_02.T), Ad_3))

Xa = np.array(Xa[:][:]).squeeze(axis=2)
Xb = np.array(Xb[:][:])
Xc = np.array(Xc[:][:])
Xd = np.array(Xd[:][:])

(964, 3)

Step 9: Camera positioning in 3D (four solutions)

In this step you will compute location of each cameras' optical center and its (calibrated) image center as points in 3D (world coordinate system). The next step 10 visulaizes the computed cameras' optical centers in red and image centers in green.

# camera's optical centers (for pair of cameras) as points in 3D world coordinate system
# 2x3 matrices: two rows with 3D point coordinates for the first and second camera
Ca = np.zeros((2,3))
Ca[1] = T1.T

Cb = np.zeros((2,3))
Cb[1] = T2.T

Cc = np.zeros((2,3))
Cc[1] = T1.T

Cd = np.zeros((2,3))
Cd[1] = T2.T

# calibrated/normalized image centers (for pair of cameras) as points in 3D world coordinate system
# 2x3 matrices: two rows with 3D po### Step 9: Cameras positions/orientations in 3D  (four solutions)int coordinates for the first and second camera
Qa = np.zeros((2,3))
Qa[0] = np.array([0, 0, 1])
temp = R1 @ np.reshape(a=[0, 0, 1], reshape=(3, 1)) + T1.reshape(3, 1)
Qa[1] = temp.T

Qb = np.zeros((2,3))
Qb[0] = np.array([0, 0, 1])
temp = R1 @ np.reshape(a=[0, 0, 1], reshape=(3, 1)) + T2.reshape(3, 1)
Qb[1] = temp.T

Qc = np.zeros((2,3))
Qc[0] = np.array([0, 0, 1])
temp = R2 @ np.reshape(a=[0, 0, 1], reshape=(3, 1)) + T1.reshape(3, 1)
Qc[1] = temp.T

Qd = np.zeros((2,3))
Qd[0] = np.array([0, 0, 1])
temp = R2 @ np.reshape(a=[0, 0, 1], reshape=(3, 1)) + T2.reshape(3, 1)
Qd[1] = temp.T

print (Cc)
print (Qc)

[[ 0.          0.          0.        ]
 [-0.02986823  0.02785665  0.9991656 ]]
[[ 0.00000000e+00  0.00000000e+00  1.00000000e+00]
 [-2.98682300e-02  2.78566527e-02 -8.34400242e-04]]

Step 10 $\text{(fully implemented)}$: 3D visualization of cameras and triangulated points (four solutions)

# visualization part
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(10,figsize = (10, 10))

ax10_1 = plt.subplot(221, projection='3d')
plt.title('Solution a')
ax10_1.scatter(Xa[:,0],Xa[:,1],Xa[:,2], c='b', marker='p')
ax10_1.scatter(Ca[:,0],Ca[:,1],Ca[:,2], c='r', marker='p')
ax10_1.scatter(Qa[:,0],Qa[:,1],Qa[:,2], c='g', marker='p')

ax10_2 = plt.subplot(222, projection='3d')
plt.title('Solution b')
ax10_2.scatter(Xb[:,0],Xb[:,1],Xb[:,2], c='b', marker='p')
ax10_2.scatter(Cb[:,0],Cb[:,1],Cb[:,2], c='r', marker='p')
ax10_2.scatter(Qb[:,0],Qb[:,1],Qb[:,2], c='g', marker='p')

ax10_3 = plt.subplot(223, projection='3d')
plt.title('Solution c')
ax10_3.scatter(Xc[:,0],Xc[:,1],Xc[:,2], c='b', marker='p')
ax10_3.scatter(Cc[:,0],Cc[:,1],Cc[:,2], c='r', marker='p')
ax10_3.scatter(Qc[:,0],Qc[:,1],Qc[:,2], c='g', marker='p')

ax10_4 = plt.subplot(224, projection='3d')
plt.title('Solution d')
ax10_4.scatter(Xd[:,0],Xd[:,1],Xd[:,2], c='b', marker='p')
ax10_4.scatter(Cd[:,0],Cd[:,1],Cd[:,2], c='r', marker='p')
ax10_4.scatter(Qd[:,0],Qd[:,1],Qd[:,2], c='g', marker='p')

plt.show()

Step 11: Reprojection errors

# Randomly select N=50 matches (pairs of features in two images) from the set of inliers for E
N = 50
ind_sample2 = np.random.choice(num_inliers, N, replace = False)

# Indicate (E) inlier matches in image 1 and image 2
plt.figure(11,figsize = (10, 4))
ax11_1 = plt.subplot(121)
plt.imshow(imL)
plt.plot(ptsL[ind[E_inliers][ind_sample2], 0], ptsL[ind[E_inliers][ind_sample2], 1], 'ob')
ax11_2 = plt.subplot(122)
plt.imshow(imR)
plt.plot(ptsR[ind[E_inliers][ind_sample2], 0], ptsR[ind[E_inliers][ind_sample2], 1], 'ob')

# project reconstructed 3D points onto both images and display them in red color
# a. convert correct points (Xa, Xb, Xc, or Xd) to homogeneous 4 vectors
# b. project homogeneous 3D points (onto uncalibrated cameras) using correct Projection matrices (KPw and, e.g. KPa)
# c. convert to regular (inhomogeneous) point
Xa_h = np.column_stack((Xa, np.ones(shape=(Xa.shape[0] ,1))))
Xb_h = np.column_stack((Xb, np.ones(shape=(Xb.shape[0] ,1))))
Xc_h = np.column_stack((Xc, np.ones(shape=(Xc.shape[0] ,1))))
Xd_h = np.column_stack((Xd, np.ones(shape=(Xd.shape[0] ,1))))

Xa_l = np.matmul(np.matmul(K, Pw), Xa_h.T).T
Xa_r = np.matmul(np.matmul(K, Pa), Xa_h.T).T

Xb_l = np.matmul(np.matmul(K, Pw), Xb_h.T).T
Xb_r = np.matmul(np.matmul(K, Pb), Xb_h.T).T

Xc_l = np.matmul(np.matmul(K, Pw), Xc_h.T).T
Xc_r = np.matmul(np.matmul(K, Pc), Xc_h.T).T

Xd_l = np.matmul(np.matmul(K, Pw), Xd_h.T).T
Xd_r = np.matmul(np.matmul(K, Pd), Xd_h.T).T


ptsL_proj = np.zeros((N,2))
ptsR_proj = np.zeros((N,2))


ax11_1.plot(ptsL_proj[:,0], ptsL_proj[:,1], '.r')
ax11_2.plot(ptsR_proj[:,0], ptsR_proj[:,1], '.r')

ax11_1.plot(Xa_l[:,0], Xa_l[:,1], '.r')
ax11_1.plot(Xb_l[:,0], Xb_l[:,1], '.g')
ax11_1.plot(Xc_l[:,0], Xc_l[:,1], '.m')
ax11_1.plot(Xd_l[:,0], Xd_l[:,1], '.k')
ax11_2.plot(Xa_r[:,0], Xa_r[:,1], '.r')
ax11_2.plot(Xb_r[:,0], Xb_r[:,1], '.g')
ax11_2.plot(Xc_r[:,0], Xc_r[:,1], '.m')
ax11_2.plot(Xd_r[:,0], Xd_r[:,1], '.k')

plt.show()

Question: how different are projected points for $SfM$ solutions a, b, c, and d? Explain.

Answer:

They are signifcantly different due to translation and rotation of the points.